﻿/*
   Created 2009 March 23 Randall Maas
   Copyright (c) 2009, Randall Maas
   All rights reserved.

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions are met:

   Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.

   Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

   Neither the name of Blackwood Designs, LLC nor the names of its contributors
   may be used to endorse or promote products derived from this software without
   specific prior written permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
   THE POSSIBILITY OF SUCH DAMAGE.
*/

/*
   This implements David Robinson's equal loudness algorithm for estimating
   the dB of sound *as heard by the ear*
    http://replaygain.hydrogenaudio.org/calculating_rg.html
    http://replaygain.hydrogenaudio.org/equal_loudness.html
 */

using System ;
using System . Text;


partial class GainControl
{
   /// <summary>
   /// The length of the dB window (seconds)
   /// </summary>
   const double dBWindow = 0.750;

   #region The current state
   /// <summary>
   /// Rolling window of the last 3/4 second of sound
   /// </summary>
   double[] MSQueue;

   /// <summary>
   /// The index into MSQueue
   /// </summary>
   int QIdx = 0;

   /// <summary>
   /// The sorted version of MSQueue
   /// </summary>
   double[] SortedQ;

   /// <summary>
   /// The index in SortedQ that nominally corresponds to the 95%
   /// </summary>
   int Q95Idx;

   /// <summary>
   /// The equal-loudness yulewalk coefficients
   /// </summary>
   double[] YuleCoefs;

   /// <summary>
   /// The equal-loudnes 150HP filter coefficients
   /// </summary>
   double[] HPCoefs;

   // The yulewalk filter state
   // Could use a buffer of 10 values and a modulus to circulate around them
   // but that is more work than just tracking them all for each sample window
   double[] YuleTmp;

   // The analysis's butterworth filter state variables
   double GA_O_1=0.0, GA_O_2=0.0, GA_I_1=0.0, GA_I_2=0.0;
   #endregion

   /// <summary>
   /// Analyzes the current sample window and estimates it dB.
   /// </summary>
   /// <remarks>
   /// This dB estimate is made using using an algorithm similar to the presented at:
   /// http://replaygain.hydrogenaudio.org/calculating_rg.html
   /// 
   /// First, we normalized the sample to approximate what the human ear perceives,
   /// using a filter described at:
   /// http://replaygain.hydrogenaudio.org/equal_loudness.html
   /// 
   /// Next, the Mean-Square of sample window is computed and stored.
   /// The Mean-Square based on the past 750 milliseconds (as defined by 
   /// dBWindow) is retrieved.  (If the current sample window is much softer,
   /// its Mean-Square value will be used instead)
   /// 
   /// Finally, an approximate dB value is computed from the Mean-Square at the 95%
   /// of the samples for the last 750ms.  If the current sample is especially quiet,
   /// the samples mean-square value overrides.
   /// 
   /// Note: David Robinson -- at the above links -- calibrates his stereo system
   /// to have computed dB provide a predictable, meaningful value.  The dB
   /// here is not calibrated...)
   /// </remarks>
   /// <param name="Samples">The buffer of sound samples</param>
   /// <param name="Ofs">The start of the samples in the buffer (before this are old samples)</param>
   /// <param name="MaxGain">The maximum gain that can be used without clipping</param>
   /// <returns>The estimated dB of the current sound (not a calibrated SPL db)</returns>
   internal double Analyze(Int16[] Samples, int Ofs, out double MaxGain)
   {
      int NumSamples = Samples . Length - Ofs;
      int MaxSample = 0;
      double Sum = 0.0;

      #region Normalize samples to Human Perception (Filter Block)
      for (int L = Samples.Length, N = Ofs; N < L; N++)
      {
         int    _V= Samples[N];
         double V = _V;
         if ( _V > MaxSample) MaxSample=  _V;
         if (-_V > MaxSample) MaxSample= -_V;

         #region Yulewalk filter
         // Perform yulewalk filter
         double S = V * YuleCoefs[0];
         for (int J=N-1,I = 1; I < 11; I++, J--)
          S += Samples[J] * YuleCoefs[I] ;
         for (int J=N-1,I = 11; I < 21; I++, J--)
          S -= YuleTmp[J] * YuleCoefs[I] ;
         #endregion

         // Store for the feedback into the next stage of the yule walk
         YuleTmp[N] = S; 


         #region High-pass filter
         // Perform butterworth filter stage, using S as an input
         double Accum = 
            (S+GA_I_2) * HPCoefs[0]
            +  GA_I_1  * HPCoefs[1];
         GA_I_2 = GA_I_1;
         GA_I_1 = S;
         Accum = Accum - GA_O_1 * HPCoefs[2] - GA_O_2 * HPCoefs[3];

         GA_O_2 = GA_O_1;
         GA_O_1 = Accum;
         #endregion


         // The square of the filtered results
         Sum += Accum * Accum;
      }

      // Copy the intermediate yulewalk state for the next
      // (this is needed since we are looking at a fairly small time window)
      for (int I = 0, J = YuleTmp.Length-10; I < 10; I++)
       YuleTmp[I] = YuleTmp[J];
      #endregion

      if (0 == MaxSample)
        MaxGain = 0;
       else
        MaxGain = 32760.0/MaxSample;

      #region Statistical Analysis to find typical dB
      // The mean square of the filtered results
      double MS = Sum / NumSamples;

      // Store the sound level the array
      MSQueue[QIdx++] = MS;
      if (QIdx >= MSQueue . Length)
        QIdx = 0;

      // Sort the buffer so we can find the 95% percentile
      Array . Copy(MSQueue, SortedQ, MSQueue.Length);
      Array . Sort(SortedQ);

      // Return the 95% 
      double X = SortedQ[Q95Idx];

      for (int I = Q95Idx +1; X < 400.0 && I < SortedQ . Length; I++)
       X = SortedQ[I];

      if (MS < X * 0.40 && MS < 12800.0)
        X = MS;

      // Convert the Mean-Square to dB by taking its root and logging it
      // The epsilon prevents log of 0
      return 10.0 * Math . Log10 (X * 0.5 + double . Epsilon);
      #endregion
   }
}
